Gravity-driven Dense Granular Flows 
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We report and analyze the results of numerical studies of dense granular flows in two and three 
dimensions, using both linear damped springs and Hertzian force laws between particles. Chute 
flow generically produces a constant density profile that satisfies scaling relations suggestive of a 
Bagnold grain inertia regime. The type of force law has little impact on the behavior of the system. 
Bulk and surface flows differ in their failure criteria and flow rheology, as evidenced by the change 
in principal stress directions near the surface. Surface-only flows are not observed in this geometry. 
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Understanding the behavior of granular materials has 
been a great challenge to scientists [1,2] and engineers 
[3,4]. One major hurdle has been the lack of a formal 
connection between the complicated but relatively well- 
understood world of contact mechanics [5], which de- 
scribes the nature and dynamics of intergranular interac- 
tions, and empirical continuum models that describe the 
macroscopic behavior of the system. 

We aim to isolate the essential features of granular 
flow, unencumbered by complicated boundary effects. To 
this end, we perform simulations of granular dynamics in 
a simple geometry: gravity driven dense flow down an in- 
clined plane, denoted henceforth as "chute flow" . Several 
remarkable features emerge: 

(1) In steady-state, the packing fraction 4> remains con- 
stant as a function of depth, beyond a dilatant surface 
region a few layers thick (Fig. 1.) The compacting in- 
fluence of increasing stress due to the weight of grains 
overhead is balanced by increasing velocity fluctuations 
towards the bottom of the assembly. 

(2) Unlike Couette flows, the entire assembly is in mo- 
tion and surface-only flows are not observed. 

(3) Components of the stress tensor and the square 
of the strain rate grow linearly with depth, indicative of 
Bagnold grain-inertia behavior [6] . 

(4) Normal stresses differ from each other [7] in a sys- 
tematic way (Fig. 2), which we do not fully understand. 

We report results of large scale molecular dynamics 
simulations of chute flow in two and three dimensions 
(2D and 3D), with interparticle interactions betweeen the 
(monodisperse) spheres modeled using both damped lin- 
ear springs and Mindlin-Hertz contact forces, with static 
friction. Detailed results of the simulations will be pre- 
sented elsewhere [8]. The main obstacle for experiments 
and simulations so far has been the difficulty of reaching 
and maintaining steady state. Previous simulations [9] 
employed very few particles or did not reach steady state 
[10]. All of these simulations were in 2D with the excep- 
tion of Walton [9]. Experiments on chute flow [11,12] did 
not involve deep assemblies. Different effects of flow were 



also studied in simulation, such as size segregation [13]. 




150 



100 



50. 75. 100. 

Distance from bottom z/d 



0.6 



2 0.5 



0.4 



0.3 



40. 



20. 



1 1 


' 9=20° ' 




^ o. 






9=26° 


tv 


1 1 




9=26° 
_24° 






_22° 
20° 



10. 



20. 30. 
Distance from bottom z/d 



40. 



50. 



FIG. 



1. Density and velocity profiles as a function of height 
from the bottom of the assemblies for 2D (top) and 3D (bot- 
tom) simulations. Inset shows a schematic of the geometry. 
Results are for 10 000 particles in 2D and 8 000 in 3D. The 
characteristic velocity un - 
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FIG. 2. Tilt dependence of the packing fraction (top) and 
the normal stress anomaly (a xx — a zz )/a zz (bottom) below 
the transitional surface layer. The dashed lines are fits to the 
forms Eq. (8 -9). The solid lines are guides to the eye. 

The 3D simulation cells contain spheres of diameter d 
and mass m, supported by a fixed bottom on the x — y 
plane. The bottom wall is constructed from a cross- 
section of a random close packing of identical spheres, 
providing a rough surface. Periodic boundary conditions 
are imposed along x and y directions. 2D simulations 
follow the same procedure, except that particles are re- 
stricted to the x — z plane and the bottom consists of a 
regular array of particles of diameter 2d. In both cases, 
there is no slip observed at the bottom. (There is some 
slip in 2D with a regular array of particles of diameter d 
at the bottom.) The gravity vector g is rotated by the tilt 
angle 6 away from the (—z) direction in the x — z plane, 
so that the free surface is normal to the z axis. In 3D, 
most of our results are for 8 000 particle systems with a 
simulation cell of size L x = 18. 6d and L y = 9.3d, result- 
ing in an assembly roughly 40 particles deep at rest. In 
2D, L x — 100g? and the number of particles varied from a 
few hundred to 20 000. In this Letter, we present results 
for N = 10 000, with a depth of about 100 particles. 

We use the contact force model of Cundall and Strack 
[14]. Static friction is implemented by keeping track of 
the elastic shear displacement throughout the lifetime of 
a contact. For two contacting particles at positions ri 
and r 2 , with velocities v 12 and angular velocities o>i, 2 , 
the force on particle 1 is computed as follows: The nor- 
mal compression <5, normal velocity v„, relative surface 
velocity v t , and the rate of change of the elastic tan- 
gential displacement u t , set to at the initiation of a 
contact, are given by 

5 = d-\\r 12 \\, (1) 

V„ = (V12 • f 12 )f 12, (2) 

v t = V12 - v„ - (wi + w 2 ) x 1-12/2, (3) 



where r i2 = ri -r 2 , fi 2 = ri 2 /||ri 2 ||, and v i2 = vi - v 2 . 
The second term in Eq.(4) arises from the rigid body 
rotation around the contact point and assures that u t al- 
ways remains in the local tangent plane of contact. Nor- 
mal and tangential forces acting on particle 1 are given 

by 



F„ = k n f(5/d) (5f\2 - t„v„) 
Ft - k s f(S/d) (-u t - r.vt) . 



(5) 
(6) 



where k n>s and r„. s are elastic constants and viscoelastic 
relaxation times respectively; f(x) = 1 for damped lin- 
ear springs or f(x) = ^fx for Hertzian contacts between 
spheres. A local Coulomb yield criterion, F t < nF n , 
is satisfied by truncating the magnitude of u t as neces- 
sary. Thus, the contact surfaces are treated as "stuck" 
while F t < fJ,F n , and slipping while the yield criterion 
is satisfied. This "proportional loading" approximation 
[15] is a simplification of the much more complicated and 
hysteretic behavior of real contacts [16]. The force on 
particle 2 is determined from Newton's third law. Each 
particle is also subject to a body force 



Fbody = mg(—zcos6 + xsin9). 



(7) 



All results are given in terms of non-dimcnsionalizcd 
quantities: Distances, times, velocities, forces, elastic 
constants and stresses are reported in units of d, tn = 
\fdfg~, vo = \fgd, Fo = mg, fen = mg/d and <7rj = mg/d 2 , 
respectively. The summary of parameters used in the 
simulations are shown in Table I. In these units, the cor- 
rect elastic constant for glass spheres with d = lOOfim 
would have been fc| lass /fcn ~ 3 x 10 10 , which would have 
prohibited a large-scale simulation. We use k n /k — 
2 x 10 5 , while controlling the coefficient of restitution for 
binary collisions through t„, assuming that we remain 
sufficiently close to the k n — > 00 limit of small deforma- 
tions. Simulations for k n /ko = 2 x 10 4 and 2 x 10 6 in 2D 
gave essentially the same results, supporting this assump- 
tion. For Hertzian contacts, the ratio k s /k n depends on 
the Poisson ratio of the material [16], and is about 2/3 
for most materials. We use the value k s /k n = 2/7, since 
this makes the period of normal and shear contact oscil- 
lations equal to each other in the damped linear springs 
case [17]. The collisional dynamics are not very sensitive 
to the precise value of this ratio. 
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TABLE I. Simulation parameters for results shown. The 
normal coefficient of restitution for 2D simulations is 0.92; it 
is a function of collision velocity for Hertzian contacts in 3D 
simulations. 
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The equations of motion for the translational and ro- 
tational degrees of freedom were integrated with either 
a third-order Gear predictor-corrector or velocity- Verlet 
scheme with a time step St = 1 x 1CP 4 . Typically 
it was necessary to run between 5 and 20 x 10 6 St to 
reach steady state, particularly when starting from a non- 
flowing state. All coarse grained quantities have been 
averaged both temporally (typically 2 to 8 x 10 6 St) and 
spatially over slices of constant z. 

The main characteristics of all the flows are: (i) The 
existence of an "angle of repose" r , such that granular 
flows can not be sustained for < r , (ii) a steady-state 
flow with a packing fraction independent of depth for 
8 r < 8 < 6 max , and (iii) for 9 > f? max , development of a 
shear thinning layer at the bottom of the assembly that 
results in lift-off and unstable acceleration of the entire 
assembly. For very thin assemblies, less than about 20 
layers, the value of 9 r depends on their depth, in agree- 
ment with experiment [12]. Here we consider only deep 
assemblies where the value of r is independent of depth, 
and focus our attention on region (ii). As seen in Fig. 1, 
the packing fraction <fi remains constant as a function of 
depth, away from the free surface and the bottom wall. 
Its value is shown as a function of 9 in Fig. 2. Results in 
2D for systems of size 5 000 and 20 000 demonstrate that 
the thicknesses of the boundary layers at the bottom and 
top are independent of the height of the assembly The 
data suggest an approximate tilt dependence for <f> of the 
form 

<p2D(e)^<p^-c 2D (0-e rt 2D) 2 , (8) 

&Zj(0)«$»S X -C3Zj(0-0r,3Zj). (9) 

where $fg x = 0.810(5) and = 0.595(5). 

In 2D, upon lowering the tilt angle below 9 r , we ob- 
serve a compaction to a polycrystalline triangular lattice 
with <p2D ~ 0.9. This causes considerable hysteresis in 
2D simulations as 9 is subsequently increased beyond 9 r : 
There is no flow until a maximum angle of stability is 
exceeded. Initial failure always occurs at the bottom, fol- 
lowed by movement of a dilation front towards the top. 
Once the system reaches steady state, 9 can be reduced 
and the system continues to flow while 9 > 9 r . On the 
other hand, in 3D there is no jump in <j> as the system 
comes to a stop, and no detectable hysteresis as the sys- 
tem is stopped and restarted by varying 9. Thus, Mohr- 
Coulomb analysis of the stress tensor [3] can be used to 
relate the flow rheology near 9 r to the Coulomb failure 
criterion as follows: The Mohr circle is the set of normal 
and shear stresses (a and r, see inset in Fig. 3) associ- 
ated with all possible shear planes. The points A and 
B at coordinates (cr zz , a xz ) and (o xx , — cr xz ) in the (<r, r) 
plane form a diameter of the circle. At a given tilt angle, 
o zz and a xz are determined by force balance, which pins 
down the location of point A (<r xz /a zz = tan6>; see Fig. 3, 
inset.) 
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FIG. 3. Left: The gap angle 89 (defined pictorially in 
the inset) as a function of height for 6 = 20° (circles), 22° 
(squares), 24° (diamonds), and 26° (triangles). The lines are 
fits that decay exponentially from S9 Buri to S9 hulii with a typi- 
cal decay length of 1.5 to 2.5d, inidcating a surface layer about 
5d to 8d thick. Right: Although <56» surf vanishes at 9 = 9 r , 
S9 hulk remains finite, suggesting that the initiation of flow is 
controlled by the surface rather than the bulk (see text.) 

However, a xx , thus the location of point B, remains in- 
determinate from these considerations. The "gap angle" 
S9 = 9mc — 9 is a measure of the difference between r/er 
along the slip plane parallel to the surface (= tan 9) and 
the largest value of rja experienced by the system (equal 
to the slope tan 9mc oi the line that passes through the 
origin and that is tangent to the Mohr circle.) For an 
ideal Coulomb material with a uniform yield criterion, 
89 = at 9 = r , when the static system is at incipient 
yield. The behavior of 80 as a function of depth at differ- 
ent tilt angles, as shown in Figure 3, reveal that: (i) SO 
becomes independent of depth below a transitional sur- 
face layer about 5ri to 8ri in thickness, (ii) the gap angle 
remains finite at the bottom and in the bulk but vanishes 
at the top surface when = r . It thus seems that al- 
though the bulk of the system has the intrinsic capability 
to withstand slightly larger tilt angles, the destabilization 
of the surface at 6 = 9 r is ultimately responsible for the 
failure and initiation of flow in the entire system. Note 
that the transitional surface layer is not directly related 
to the dilatant layer; it is much thicker near = r and 
penetrates well into the region of constant density. In- 
terestingly, in 2D the gap angle at r is actually larger at 
the surface compared to the bulk, and both gap angles 
remain finite. However, the presence of hysteresis pre- 
cludes us from applying the Mohr-Coulomb analysis in 
2D. 

A feature that distinguishes granular flows from New- 
tonian fluid flows is that normal stresses, i.e., diagonal 
terms in the stress tensor, are in general not equal to 
each other [7]. Although a xx w a zz in our simulations, 
we observe small but systematic deviations, which are 
depicted in Fig. 2 by plotting the normal stress anomaly 
in the bulk, (a xx — <j zz )/<j zz , as a function of tilt an- 
gle. These deviations are likely to be due to a constitu- 
tive equation of the flow rheology which we have not yet 
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been able to determine. In addition, in all 3D runs, a yy 
is smaller than the other normal stresses by 10 — 15%, 
suggesting that consolidation and compaction normal to 
the shear plane is poorer. 

Another question of particular interest is the relation- 
ship between the stresses and strain rate. Shear stress a xz 
is a linear function of strain rate 7 = dv x /dz for viscous 
fluids, and quadratic in 7 for granular systems in the Bag- 
nold grain-inertia regime. The latter result is rather gen- 
eral: When typical stresses on grains are large compared 
to the weight of individual particles but not large enough 
to significantly distort the spheres (1 <C ct/cto <C k n /ko), 
the only relevant time scale is 7 , which forces a xz oc j 2 
simply by dimensional analysis. As an example, Fig. 4 
shows the relationship between shear strain rate 7 and 
shear stress a xz for 2D and 3D cases. Below the transi- 
tional surface layer, and away from the bottom wall, both 
systems exhibit Bagnold scaling, indicated by the dashed 
lines in Fig. 4. Data for 2D suggest an offset of about 
1.5<7o in the stress, possibly due to corrections from the 
body force on individual spheres. Such an offset is not 
needed for an acceptable fit to the 3D data. 
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FIG. 4. Strain rates plotted as a function of shear stress a xz 
for 2D (top) and 3D (bottom) reveal Bagnold scaling (dashed 
lines) in the constant packing fraction regime away from the 
free surface and the bottom wall. For 2D, an offset of 1.5cto 
in the shear stress was needed for an acceptable fit. 

In order to probe the sensitivity of the results to the 
exact form of the stiff elastic response, we have also 
performed runs in the 3D system with linear damped 
springs, keeping all other parameters fixed [8]. Remark- 
ably, the packing fraction profiles and the normal stress 
anomaly remained virtually the same, and strain rate 
profiles were changed only by a global factor of about 



1.35, suggesting that results are not too sensitive to the 
particular force scheme selected. 

The lack of a regime with only surface flow is in con- 
trast with experimental observations of avalanche flows 
in rotating drums [2] . Since periodic boundary conditions 
are used, our simulations correspond to an infinitely large 
system with finite depth and fixed surface tilt, and there- 
fore constitutes a different system. Experiments on chute 
flows [7] also lack a regime of surface only flow; although 
this might be attributed to side- wall friction that necessi- 
tates higher tilt angles to initiate flow. Although there is 
no fundamental reason we can find that prohibits surface- 
only flows in this geometry, it appears that they are un- 
likely to occur in an assembly of monodisperse spheres. 
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